function res = cctangent(c1, c2)
    %c1 and c2 are the centers of gears 1 and 2 in the form [c1x, c1y, r1]
    hold on
    r1 = c1(3);
    r2 = c2(3);
    
%    a = sqrt( (c2(2) - c1(2))^2 + (c2(1) - c1(1))^2 );
%    theta = asin(r2/a);
    m =  ( c2(2) - c1(2) ) / ( c2(1) - c1(1) )
    x1 = r1 * sqrt( 1 / ( 1 + (1/-m)^2 ));
    y1 = 1/(-m) * x1;
    x2 = r2 * sqrt( 1 / ( 1 + (1/-m)^2 ));
    y2 = 1/(-m) * x2;
    t1 = [c1(1) + x1 , c1(2) + y1];
    t2 = [c2(1) + x2 , c2(2) + y2];

    %p1 is the tangent point on the circle centered on c1 with r = r1 - r2
    res = [t1, t2];
    %res = cctangeqr(p1, c2, r2);
    
    
    line([t1(1),t2(1)],[t1(2), t2(2)]);
    plot(c1(1),c1(2),'ro');
    circle(c1, r1,1000);         %Large circle around c1
    circle(c2, r2,1000);         %Large circle around c2
    circle(c1, r1-r2, 1000, 'r');
end



function res = cctangeqr(p1, p2, r)
    plot(p1(1),p1(2),'ro')
    plot(p2(1),p2(2),'ro')
    
    m = ( p2(2) - p1(2) ) / ( p2(1) - p1(1) );
    x = r * sqrt( 1 / ( 1 + (1/-m)^2 ));
    y = 1/(-m) * x;
    t1 = [p1(1) + x, p1(2) + y];
    t2 = [p2(1) + x, p2(2) + y];
    res = [t1, t2];
    
    plot(t1(1),t1(2),'ro')
    plot(t2(1),t2(2),'ro')
    circle(p1,r,1000)           %small circle around c1
    circle(p2,r,1000)
    line([t1(1),t2(1)],[t1(2), t2(2)])
end